Finite-Difference Investigation of Axisymmetric 
Inviscid Separated Flows with Infinitely-Long 
Cusp-Ended Stagnation Zone. Flow around a Sphere 



M. D. Todorov 

Dept. of Differential Equations, Institute of Mathematics and Informatics, 
Technical University of Sofia, Sofia 1756, Bulgaria 
e-mail: mtod@vmei.acad.bg 



Abstract 

The classical Helmholtz problem is applied for modelling the axisymmetric in- 
viscid cusp-ended separated flow around a sphere. Two coordinate systems are 
employed: polar for initial calculations and parabolic the latter being more suit- 
able for investigation of infinitely long stagnation zones. Scaled coordinates are 
introduced and difference schemes for the free-stream equation and the Bernoulli 
integral are devised. The separation point is not initially prescribed and is defined 
iteratively. A separated flow with vanishing drag coefficient is obtained. 



1. Introduction 

In an attempt to explain the existence of a sizable drag force upon a submerged 
body even for vanishing viscosity, Helmholtz [IIJ introduced the notion of discontinuous 



ideal flow consisting of a potential and stagnant parts; these matching at an unknown 
stream surface. The idea of discontinuous ideal flow was successfully applied by Kirchhoff 
[TT|j for bodies with sharp edges and later developed by Levi-Civita ||12|| , Villat [^0 



Brodetsky [|], etc. for bodies with curved profile when additional condition for smooth 
separation (Brillouin- Villat condition) is to be satisfied. All these solutions are planar 
and based on the hodograph method. Unfortunately this powerful tool is not capable 
for ideal flows characterized by axial symmetry. Therefore the efforts in solving of such 
kind flows is mainly confined to the numerical approach. Hitherto there are known several 
approximate methods for study of axisymmetric ideal flows. The most important methods 
appear to be: the integral one used at first by Trefftz |19| and later extended by Struck 



16fl ; the relaxation one applied by Southwell&Vaisey |1| and developed by Brennen |§ 



(for detailed reference see |2T|, ||). Similarly to the plane flows in the case of curve bodies a 
smooth separation condition or any else semi-empirical assumptions are suggested in order 
to yield satisfactory forecast concerning the velocity and pressure distribution, detachment 
point and drag coefficient [|, 13|. Particularly Southwell&Vaisey by working in the 



physical plane obtained only cusp-ended cavity behind a sphere. We also calculated such 



kind stagnation zone behind a sphere [H] by means of finite-difference scheme at that 
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without pre-conditioning the separation point. Now we aim at utilizing the improved 
difference scheme, which was developed and applied for the planar inviscid flow around 
circular cylinder |18j for investigation of a separated axisymmetric flow around a sphere. 



Following our approach we use two different coordinate systems: a polar spherical 
coordinate system for initial calculations and a parabolic coordinate system the latter 
being topologically more suited for solving the free-stream equation outside infinitely- 
long stagnation zones. We switch from polar coordinates to parabolic ones after the 
stagnation zone has fairly well developed and has become long enough. 

2. Posing the Problem 

Consider the steady inviscid flow past a circle - an arbitrary meridian cross section of 
a sphere. The direction of the flow coincides with the line 9 = 0, tt of the polar coordinates 
and the leading stagnation point of the flow is situated in the point 9 = ir. The axially 
symmetry enables to study the flow in the meridian halfplane only. 

Dimensionless variables are introduced as follows 



,/ ^ I r P-Pc fj I IT I Poo-Pc 

Y = T0TT i t — — , q = -j— — — , cr = v La , r = v Lr 



(2.1) 



where L is the characteristic length of the body (2a for a sphere of radius a) , - velocity 
of the undisturbed flow; p c - the pressure inside the stagnation zone; - the pressure 
at infinity, r - the polar radius, a, r-the parabolic coordinates, n - the cavitation number, 
which for flows with stagnation zones is equal to zero. Without fear of confusion the 
primes will be omitted henceforth. 

2.1. Coordinate Systems 

In terms of the two coordinate systems (polar spherical and parabolic) equation for 
the stream function ip reads: 

1 ( I \ 1 f ^ \ n 

:W0r + -r -— ? =0, or 




sin# r 2 ysin^ 

The undisturbed uniform flow at infinity is given by 

r 2 ^oo sin 2 9 

^Uoc ~ 2 ' ° r ^1™,™ ~ vtUoc • (2.3) 

On the combined surface "body+stagnation zone" hold two conditions. The first 
condition secures that the said boundary is a streamline (say of number "zero" ) 

i>(R(9),9) = 0, 9 G [0,tt] or i/j(S(t),t) = 0, r G (0, oo) , (2.4) 

where R(9), S(t) are the shape functions of the total boundary in spherical or parabolic 
coordinates, respectively. As usually we use the notation Fi for the portion of boundary 
representing the rigid body (the sphere) and T 2 - for the free streamline (Fig.l). 

On T 2 the shape function R(9) is unknown and it is to be implicitly identified from 
Bernoulli integral with the pressure equal to a constant (say, p c ) which is the second 
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condition holding on the free boundary. For the two coordinate systems one gets the 
following equations for shape functions R(6) or S(r): 



q + 



1 4+« 



r 2 sin 2 9 r 2 
9eT 2 



1. 



or 



T = R(d) 



q + 



2 . 2 

a + t 



1. (2.5) 



<x=S(r) 



The boundary value problem 
symmetry conditions 



r G T 2 

( |2.5| ) is completed with the additional 



-7^=0, 9 = 0,7v or ^ = , r = 



(2.6) 



In spherical coordinates along with ip it is convenient to introduce new function \I/ 
-*—E. Then the dynamical condition (|2.5|a) takes the form: 



tpr 



Obviously ^r\ r=R(e) rsm 
name \1/ stream function too 



*1 

e r 2 



i . 



r=R(d) 



r=R(6) 



r=R(8) rsini 



r=R(8)' 



(2.7) 



Without confusion we will 



2.2. Scaled Variables 

Following |], 0, [| we introduce new scaled coordinates: 

i 1 = rR-\9), v = (t-S(t), 

which render the original regions to semi-infinite strips. 

If we denote £ = 9 or £ = r depending on the particular case under consideration, 



then in terms of the new coordinates (77, £), the governing equation ( |2.2|) takes the form 
Atyjr, + B{b^ - C(^)„ - + + (/^)* = , (2.8) 

where 



sin 1 



fl' 1 

i2 sin# 



e = , / = 0; 



/ j^i \ z ^/ 

t4ee?] 2 -|-^ — , .B ee sin 6* , C ee 77— , D ee 77 sin 6 1 
V R I R 



or 



1 S" 1 

b=l, d=S', e EE ■ — . /EE-; 

r] + b r r 

A ee 1 + S' 2 , B = l, C = S' , D = l. 
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Similarly to [18| we use the "relative" function ip 

[rjR(9) sinO} 2 



^(rj, t) = tp(r}, t) - (77 + S(t))t , 



which is obviously a solution to eq.( |2.8| ) and which we loosely call stream function. The 
asymptotic boundary condition then becomes 



'0 







or 



4' 



V=Voo , T=T a 



while the non-flux condition on V transforms as follows 

[R{9) sin9] 2 



or 



1> 



?7=0 







-5-(t)t. 



(2.9) 



(2.10) 



Thus eqs.([2~8|), (|2.9|), ( f2.10| ), ( |2.6| ) define a well posed boundary value problem provided 
that functions R{9) and S(t) are known. On the other hand in the portion I"2 of the 
boundary (where these functions are unknown) they can be evaluated from the Bernoulli 
integral ( |2.5| ) and ( |2.7|) which now becomes an explicit equation for the shape function 



R 2 + R' 2 
R G sin 2 9 

R 2 + R' 2 
R A 



dip 



drj 



+ R 2 {9) sm 2 9 



rj=l 



< e < e* 



+ R(8) sin6 



/2 



1 + S< 

S 2 + T 2 



dtp 
drj 



r)=0 



1, 

or 
1 , 



(2.11) 



T* < T < OO 



Here (77, 0) = *(/?, 9) - rjR{9) sin 9. 



3. Forces Exerted on the Body 



Apparently the presence of a stagnation zone breaks the symmetry of the integral 
for the normal stresses and hence D'Alembert paradox is not hold. If denote by n the 
outward normal vector to the sphere £ and by da - the surface element of the sphere, 
then the force acting upon the body is given by 



R 



<j> pnda 



(3.1) 



It is not difficult to obtain for the drag and lifting-force coefficients of every meridian 
cross section the following expressions 

C x = - I\r{9) sin(0) \R{9) cos 9 + R'{9) sin 9} d9 or C x = f qS(r)r [S{t) + S'{t)t] dr 
J,r J " (3.2) 



C I, 



0, 
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where the dimensionless pressure is given by 



q = 1 
or 

Q = 1 



R 2 + R* 
i? 4 

1 + S' 2 
S 2 + r 2 



dr] 



+ R sin 9 



rj=l 



dip 



dr) 



+ T 



ri=0 



(3.3) 



4. Difference Scheme and Algorithm 



4-1. Splitting scheme for the free-stream equation 

The computational domain being infinite is reduced to finite one after appropriately choos- 
ing the "actual infinities" . In order to take into consideration the topological and dynamic 
features of the flow we employ non-uniform mesh, which was presented in detail at 



Let us denote the spacings of the mesh by h i+ i 



Vi+l - Vi ) 1 



1,... ,M and 



gj + i = — £j , j = 1, . . . , N. We solve the boundary value problem iteratively using 
the method of splitting of operator. Upon introducing fictitious time we render the 
equation to parabolic type and then employ the so-called scheme of stabilising correction 
22j. On the first half-time step we have the following differential equations (At is the 



time increment) 



= B fi A 2 (6A 2 ^ + ') . + A ij A 1 {A 1 i; n ) ij - 6^1^% 
-D ij A 2 (dAiip n ) ij + A l (eip n ) ij + A 2 (fi> n ) ij 



(4.1) 



for i = 2, • • • , M, j = 2, . . . , N 

The second half-time step consists in solving the following differential equations 



c +1 



A ii (A 1 (A 1 r +1 ) ij -A 1 (A 1 ^) i 



(4.2) 



for % = 2, . . . , M, j = 2, . . . , N. The last two equations (4.1)-(4.2) are completed with 
respective boundary conditions [0]. Here 

Ai(.) - = ^(-k + 0(^+0, 



A 2 



9 

si 1 



are the usual difference operators based on three-point patterns with second order of 
approximation. 

Thus the b.v.p. for the stream function is reduced to consequative systems with sparse 
(three-diagonal) matrices, which are solved iteratively 0. 
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Since the condition for numerical stability of the elimination is not satisfied for all 
points of domain here a "non- monotonous progonka" (see fT4| , [|) is employed like at [|18| . 



4-2. Difference Approximation for the Free Boundary 



Following [HJ in the present work we use the dynamic condition ( |2.5| ) in spherical 
coordinates only, so that we present here just the relevant scheme in spherical coordinates. 
The equations ( |2.11| ) can be resolved for the derivative R'{0) when the following conditions 
are satisfied: 



def R\6) sin 2 6 ^ d$ 
Q{9) = T 2 >1; T= ^ 



+ (R(9) sin i 



or 



Q{0) 



def R 2 {9) 



T - — 

T 2 ' Tfy 



+ R(8) sin6. 



(4.3) 



r,=l 



The above inequalities are trivially satisfied in the vicinity of the rear-end stagnation 
point inasmuch as that for 9 — > one has T — > or T — > and hence — 



T 2 



oo or 



Y2 — > oo. The first inequality, however, is indeterminated at the point 6* = due to the 



ratio 



sin 9 
T(9)- 



For the shape function Rj of free line is solved the following difference scheme 

Rj + Rj-i 



Rj-i-Rj 



9j- 



\ 



(Rf^smOj 



+ 



(Rf^ysme^ 



- 1 



or 



- 9j 



Rj + Rj-i 



\ 



(4.4) 



+ 



for j = j*,... ,2 , whose approximation is O(gj). Only in the detachment point the 
difference scheme is different, specifying in fact the initial condition, namely 

^R{e*) + R r 



Rj, - R{6*) = g* 



\ 



'(Rf^sindj 



J-j* 



(R 2 (6*) sin 0*' 



or 



R i 



R{6* 



,R(6*) + Rj 



\ 



R£ 



+ 



'R{9*) 



where R without a superscript or "hat" stands for the known boundary of rigid body. 

At last a relaxation is used for the shape-function of the free boundary at each global 
iteration a according to the formula: 

R a+1 = uRj + (l- uj)R° , 

where u> is called relaxation parameter. 



4-3. The general Consequence of the Algorithm 
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Each global iteration contains two stages. On the first stage, the difference problem 
for free-stream equation is solved iteratively either in polar spherical or in parabolic 
coordinates (depending on the development of the stagnation zone). 

The second stage of a global iteration consists in solving the difference problem for 
the free surface in polar spherical coordinates. 

Through the indetermination at the axis of symmetry we use the difference scheme 
( fOj a.) only during the first several iterations (in polar spherical coordinates). The cal- 
culation of the shape of the far weak (in parabolic coordinates) we carry out using the 
scheme ( |4.4| b). The latter appears to be more convenient and efficient because the loss 
of accuracy and 'numerical' instability in vicinity of the axis of cusp are avoided. The 
criterion for convergence of the global iterations is defined by the convergence of the shape 
function, namely 



max 

3 



po+1 pa 



ce+1 



< 1(T 4 . (4.5) 



The obtained solutions for the stream function and the shape function of the boundary 
are the values of the last iteration if;^ = ipt) +1 an d Rj = respectively. Then the 

velocity, pressure, and the forces exerted from the flow upon the body are calculated. 

5. Results and Discussion 

The numerical correctness of scheme (4.1), (4.2) is verified through usual experiments 
including a doubling the mesh knots and varying the 'actual infinity' We used differ- 
ent meshes with sizes M x N : 41x68, 81x136, 101x202, etc. Respectively, the actual 
infinity r]^ assumed in the numerical experiments the values 10, 20. The dependence 
of the numerical solution on the time increment At is also investigated and it is shown 



that the scheme of fractional steps for the stream function has a full approximation [p2|| . 
Comparing the different finite-difference realizations of the solution we choice the follow- 
ing 'optimal' values of the governing parameters: step of the fictitious time At = 0.5, 
relaxation u = 0.01 and 'actual' infinity rj^, = 10. 

In Fig.2-a are presented the obtained shapes of the stagnation zone behind the sphere 
and in the near wake for resolutions 41 x 68, 81 x 136 and 101 x 202 and value of relaxation 
parameter: uj = 0.01. 

Evidently the agreement among the calculated shapes of the free boundary near the 
body corresponding to these three meshes is very well. The logarithmic scale is used in 
Fig.2-b in order to expand the differences between the different solutions making them 
visible in the graph. As clearly it is shown the curves are indistinguishable till distance 
150 calibers and the relative error is less than 1%. The relative error between the meshes 
81 x 136 and 101 x 202 at distances more than 150 calibers does not exceed 4%. At the 
same time the relative error between the mesh 41 x 68 and the else two ones increases and 
reaches 7-8% at distance 200 calibers. Obviously that mesh is not enough fine and appears 
to be coarse for calculating the shape function at large distances behind the sphere. The 
obtained results warrant conclusion that the scheme is fully effective in solving the free 
boundary till 200 calibers. The very good comparison supports the claim that indeed a 
solution to the Helmholtz problem has been found numerically by means of the developed 
in the present work difference scheme. The calculated here dimensionless pressure q is 
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shown in Fig. 3. The agreement among the obtained pressure curves corresponding to 
different mesh resolutions is excellent. In the stagnation zone the pressure is in order of 
10~ 4 in accordance with the assumption that the unknown boundary is defined by the 
condition q = 0. The amplitude of the minimum of q is smaller than 1.25 the latter 
being the value for ideal flow without separation. This means that the stagnation zone 
influences the flow upstream. The calculated magnitude of the separation angle (measured 
with respect the rear end of the sphere) varies between 69.42° for mesh 41 x 68 and 69.7° 
for mesh 101 x 202. It is interesting to note that the calculated here drag coefficient C x 
has a magnitude between .5848 x 10 -3 — .5704 x 10~ 2 obtained for the different resolutions, 
i.e., we conclude that in order of approximation of the scheme C x = 0. Then similarly to 
the separated flow around circular cylinder we can name the obtained separation angle 
'critical' (see ]18[ |{|). Hence in the case of axisymmetric flow around sphere there also 
exists a inviscid separated flow for which the D'Alembert paradox holds. Trough the 
disscused features of the obtained Helmholtz flow we can assume it is an axisymmetric 
analogue of the Chaplygin-Kolscher flow around circular cylinder. 

6. Concluding Remarks 

The separated inviscid flow behind a sphere is treated as a flow with free surface 
- the boundary of the stagnation zone (Helmholtz problem). Scaled coordinates are 
employed rendering the computational domain into a region with fixed boundaries and 
transforming the Bernoulli integral into an explicit equation for the shape function. A new 
free-stream function is introduced and thus the numerical instability near the symmetry 
axis is avoided. Difference scheme using coordinate splitting is devised. Exhaustive set of 
numerical experiments is run and the optimal values of scheme parameters are defined. 
Results are verified on grids with different resolutions. The obtained here shape of the 
stagnation zone is of infinitely long cusp and respective separated flow has vanishing drag 
coefficient. The detachment point is not prescribed in advance and it is defined iteratively 
satisfying the mere Bernoulli integral there. 
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FIGURE CAPTIONS 

figi-gif 

Figure 1: Posing the problem 

sphnear.gif 
(a) behind the sphere 

sphfar.gif 
(b) far wake 
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Figure 2: The obtained separation lines for relaxation parameter u = 0.01 and different 
resolutions: 41 x 68; - - 81 x 136; - - 101 x 202. 

sphpres.gif 

Figure 3: The pressure distribution for relaxation parameter ui = 0.01 and different 

resolutions: 41 x 68; - - 81 x 136; - - 101 x 202; - - inviscid nonseparated 

flow. 
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